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Abstract 

Background: Contact-dependent inhibition (CD!) lias been recently revealed as an intriguing but ubiquitous 
mechanisnn for bacterial competition in which a species injects toxins into its competitors through direct physical 
contact for growth suppression. Although the molecular and genetic aspects of CDI systems are being increasingly 
explored, a quantitative and systematic picture of how CDI systems benefit population competition and hence alter 
corresponding competition outcomes is not well elucidated. 

Results: By constructing a mathematical model for a population consisting of CDI+ and CDI- species, we have 
systematically investigated the dynamics and possible outcomes of population competition. In the well-mixed case, 
we found that the two species are mutually exclusive: Competition always results in extinction for one of the two 
species, with the winner determined by the tradeoff between the competitive benefit of the CDI+ species and its 
growth disadvantage from increased metabolic burden. Initial conditions in certain circumstances can also alter the 
outcome of competition. In the spatial case, in addition to exclusive extinction, coexistence and localized patterns 
may emerge from population competition. For spatial coexistence, population diffusion is also important in 
influencing the outcome. Using a set of illustrative examples, we further showed that our results hold true when the 
competition of the population is extended from one to two dimensional space. 

Conclusions: We have revealed that the competition of a population with CDI can produce diverse patterns, 
including extinction, coexistence, and localized aggregation. The emergence, relative abundance, and characteristic 
features of these patterns are collectively determined by the competitive benefit of CDI and its growth disadvantage 
for a given rate of population diffusion. Thus, this study provides a systematic and statistical view of CDI-based 
bacterial population competition, expanding the spectrum of our knowledge about CDI systems and possibly 
facilitating new experimental tests for a deeper understanding of bacterial interactions. 

Keywords: Contact dependent inhibition. Bacterial population. Competition, Extinction and coexistence. Spatial 
aggregation 



Background 

Bacteria are highly social and present dominantly in 
the form of complex communities where they interact 
through a variety of fashions [1-3]. Among all types of 
bacterial interactions discovered, competition has been 
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identified as the most prevalent by recent studies [4-6]. 
The ubiquity of competition has been mainly attributed to 
the limited space and resources of natural environments. 
To maximize their survival and reproduction, bacteria 
have indeed developed numerous competition strate- 
gies, including interference competition where a species 
directly harms another via active production of toxin and 
other effectors [4,6-8]. 

Interference competition was initially shown to be 
mediated by diffusible soluble factors, such as antibi- 
otics and bacteriocins. These effector molecules serve to 
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potently decrease survival and reproduction of neighbor- 
ing bacteria at a long range spatial scale [2,9-11]. For 
instance, Lactococcus lactis produces and secretes nisin, 
a small antimicrobial peptide, into the extracellular mil- 
ium (e.g., milk) to efficiently inhibit other bacterial species 
for lactose competition [12]. A set of recent studies, how- 
ever, have illustrated that interference competition can 
also occur through direct physical contact between cells, 
revealing a new class of competition with an interaction 
scale restricted to nearest neighbors [13-17]. 

Furthermore, studies have uncovered a surprisingly high 
degree of diversity among these contact-dependent inhi- 
bitions (GDIs). They occur across a wide range of organ- 
isms including both Gram- negative and positive bacteria 
[17-21], use different toxins similar to nuclease, tRNAse 
and DNase, and further exploit various delivery machiner- 
ies spanning Type III, IV, V, and VI secretion systems 
[13-15,22,23]. Moreover, it has been shown that certain 
strains even have multiple GDI modules and multiple 
toxins for competition [17,24]. 

Despite their structural and compositional diversity, all 
of the GDI systems possess a common mode of action in 
which growth inhibition toxins are deployed into competi- 
tor cells via direct cell to cell physical contact. A repre- 
sentative example is the GdiBAI system discovered in the 
enterobacterium E. coli EG93 [14]: The system consists 
of three functional components: GdiA, the toxin effector, 
GdiB, the ^-barrel protein localized to the outer mem- 
brane for effector export, and Gdil, the immunity protein. 
Upon contact with a neighboring cell, a GDI equipped cell 
employs GdiB to inject toxin GdiA into the target cell to 
inhibit its growth while expressing the immunity protein 
Gdil to prevent autoinhibition. 

The intriguing functionality and characteristics of GDI 
motivate us to ask the following questions: How does 
GDI impact the competition between a GDI equipped 
(GDI+) and deficient (GDI-) species? Gan the two species 
coexist, or does extinction always occur for one of the 
species? Due to the intrinsic association of protein expres- 
sion and metabolic load, will the growth disadvantage of 
GDI+ species counteract its competition advantage from 
toxin production? How does cell motility alter the compe- 
tition outcome in different environmental settings, such 
as liquid or solid agar? 

Although current experimental efforts [13,23-26] have 
started to address some of the above questions, a system- 
atic and quantitative understanding of GDI-based bacte- 
rial competition has not been achieved. In particular, it 
is not clear how the competition outcome is influenced 
by the inhibition advantage of GDI-based competition 
and the growth disadvantage associated with metabolic 
load. Moreover, most experimental efforts have primar- 
ily focused on liquid culture settings where populations 
are well mixed [14,24,25] but little has been elucidated 



when competitions occur in space. There is hence a clear 
need for a systematic investigation of GDI-based bacterial 
competition that integrates the tradeoff between com- 
petitor inhibition and metabolic cost with spatiotemporal 
dynamics. 

Here, we present a mathematical model to describe a 
bacterial population with GDI+ and GDI- species that 
compete through both contact-dependent inhibition and 
nutrient utilization. With this model, we studied the com- 
position of the competing population in the well-mixed 
case, showing that the outcome is always extinction for 
one of the species depending on initial conditions as well 
as the tradeoff between inhibition strength and metabolic 
cost. We then continued to investigate the dynamics of the 
population in one dimensional space, revealing possible 
spatial coexistence of the species through aggregation. To 
acquire a statistical understanding of the coexistence pat- 
terns, we further conducted a systematic survey into the 
spatial structure of the competing populations by altering 
the interplay between inhibition strength and metabolic 
growth disadvantage. A set of illustrating tests in two 
dimensional space was also implemented to demonstrate 
the generality of our results. We finally conclude by 
summarizing our findings and discussing possible future 
developments. 

Results 

A mathematical model of bacterial competition with GDI 

Bacterial competition has been modeled through cou- 
pled systems of ordinary differential equations, dating 
back to the work of Lotka and Volterra [27,28]. Later 
developments of the original Lotka-Volterra model have 
incorporated the effects of nutrient limitation and 
species diffusion [29-36]. Typically, these models incor- 
porated only on-site interactions (i.e. the interaction 
range for competing species is taken to be infinites- 
imally small). For bacteria competing through GDI, 
however, it is natural to explicitly include a finite 
interaction range in order to account for the intrin- 
sic nearest-neighbor effects of toxin injection in the 
system. 

Previous studies on competition with a finite interac- 
tion range have shown that spatial aggregation is possi- 
ble in the long time limit [37-43]. For instance, a single 
species with nonlocal interactions has unstable homoge- 
neous states for certain functional forms of its growth rate, 
leading to a spatial distribution of the population with 
clumps of high species concentration separated by regions 
with low density [41,44-47]. For two species systems, it 
has been shown that systems with specific features, such 
as systems with AUee effects and completely symmetric 
interactions, may have steady spatial structures [37,40,48]. 
Despite these advances, a mathematical model appropri- 
ate for bacterial populations with the nearest neighbor and 
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asymmetric interactions implemented through GDI has 
not been investigated. 

To enable a quantitative and systematic investigation of 
CDI-based bacterial competition, we constructed a math- 
ematical model that describes a two-species population 
with one CDI+ (GDI equipped) and the other GDI- (GDI 
deficient) as follows: 

du 
dt 



u(a-u-v-c j f(\x - x\)vix)dx^ + Du'^'^u 



dv 



V) + D^V^v 



(1) 



where u and v are the concentrations of the GDI- and 
GDI+ species respectively. This model incorporates the 
asymmetrical interaction from the CDI+ to GDI- species 
through the integral term {—uc ff(\x — x\)^(x)dx)> the 
effects of limited shared resources through the logis- 
tic growth terms, and cell motilities through the corre- 
sponding diffusion terms. Here, nondimensionalization 
has been implemented with the original model detailed in 
Additional file 1: Section 1. 

Due to the intrinsic direct physical contact of GDI, 
bacterial interactions follow a discrete, nearest neighbor 
fashion, i.e.,/(|;v/ — xj\) = 1/38 when J = 1, /, / — 1 and 
fi\Xi — Xj\) = 0 otherwise. The above model can thus be 
discretized and rewritten in one spatial dimension as: 

dui 



dt 



= Ui(a - Ui - 
D 



■Cl(V/+l +V/ + V/_l)) 



dvi 



i = 1, 2, . 



- 2ui + Ui-i) 
D 



2v,- + Vi-i) 



(2) 



where 8 is the grid spacing, Ui (v/) is the concentration 
of the GDI- (GDI+) species in space /5, N is the num- 
ber of total grid points, and D is the diffusion constant 
{Du = Dy = D assumed for simplicity). In addition, we 
have used ci = c/3 to reflect the fact that there are three 
interaction locations for GDI in one dimension including 
one on-site and two nearest neighbors. Furthermore, ci 
is a constant here, corresponding to a constitutive expres- 
sion of the GDI system as seen in certain species such as 
E. coli EG93 [24]. Periodic boundary conditions are 
imposed so that u^j^i = ui and uq = um (similarly for v). 

Although this model is generic and does not describe 
biochemical details of the competing population, it cap- 
tures the main biological features of the system in realistic 
settings: (1) Gellular growth is described using a logis- 
tic growth model, which accounts for nutrient limitation 
in experimental settings and has been adopted widely by 
previous theoretical studies [35]; (2) Gellular movement 
is modeled with a random diffusion term, which is ade- 
quate to describe non-motile cells with passive diffusion 



or motile cells without chemotaxis; (3) Gonstant growth 
inhibition of the GDI system is modeled, which is sup- 
ported by experimental evidence showing that certain 
bacterial species have constitutive production of the GDI 
machinery [24]. 

Our model has four key parameters, namely a, ^, ci, 
and A which describe the growth rates of the GDI- 
and + species, inhibition strength from the GDI+ to GDI- 
species, and the motility of both species respectively. It is 
evident that if the GDI+ species (v) has a larger growth 
rate than the GDI- species (u) (i.e., a/fi < 1), the outcome 
of competition will be the extinction of the GDI- species 
due to the additive combined affects of slower growth 
and nearest neighbor inhibition. However, the outcome 
is unclear when the metabolic load of producing proteins 
involved in GDI causes a/p > 1. We are thus motivated to 
understand how the tradeoff between inhibition strength 
and growth disadvantage influences the outcome of com- 
petition and what kinds of possible competition outcomes 
may result. In the following, we consider the tradeoffs 
between relative growth advantage and inhibition over a 
range of cell motility. 

Well-mixed case 

We begin our investigation by considering the well-mixed 
case for our system, corresponding to a liquid assay or 
£) ^ 00 in the model. In this case, the solutions are 
assumed to be homogenous in space so that Ui = u and 
vi = V for all By eliminating the spatial interactions of 
the system, we get an inhibition contribution of 3c\v {cv) 
to the growth rate of u at each site. The resulting equations 
are: 



du 

—- — u(a — u — v(l + c)) 
dt 

dv 

= v(p -u-v) 

at 



(3) 



from which we find four steady states, including {u, v) = 
(0, 0), (u, v) = (a, 0), (M, V) = (0, fi), and (m, v) = {fi-{a- 
p)lc,{a-P)lc). 

The linear stability of these homogeneous steady states 
changes as we vary the parameters of the model due 
to the tradeoff between inhibition strength c and rel- 
ative growth advantage a/P (Figure lA). When the 
growth advantage of species u (GDI-), given as a/yS, 
greatly outweighs the inhibition constant c (the upper 
white region of Figure lA), species u drives v to extinc- 
tion regardless of the initial conditions as illustrated in 
Figure IG. Gonversely, when the growth advantage of 
species u is less than one (the lower white region of 
Figure lA), species v always drives species u to extinction 
(Figure ID). There does exist, however, a certain param- 
eter regime, namely 1 < al^ < 1 + c (the shaded 
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Figure 1 Well-mixed competition of a bacterial population consisting of a CDi+ and a CDi- species. (A) Phase diagram of the competing 
population. Grey region: competition of the population may lead to the extinction of one of the two species, depending on the initial conditions of 
the system. Upper white region: the CDI+ species {v) always goes extinct regardless of initial conditions. Lower white region: the CDI- species {u) 
always goes extinct regardless of initial conditions. (B) Plot of the Lyapunov function [36] of the competing population in the regime when 
extinction is possible for either species. The parameter set here corresponds to the blue cross in the grey region of Panel A. Notice that the minima 
for the system lie on the respective axes, showing two possible extinctions. The dark circle is the saddle point of the system. C-E. Sample time 
course trajectories of the two competing species. (C) Species u (CDI-) always dominates when growth advantage outweighs inhibition. Here, 
a/^ = 2 and c = 0.5, corresponding to the cross in the upper white region of Panel A. (D) Species v (CDI+) always take over the population when it 
has a growth advantage and exerts inhibition. Here, a/^ = 0.5 and c = 2.0, corresponding to the cross in the lower white region of Panel A. (E) 
Both of the species may dominate depending on their initial conditions. Here, c = 2.0 and a/^ = 2.0, corresponding to the cross in the grey region 
of Panel A. In each panel (C-E), 1 00 pairs of trajectories (orange: CDI+, blue: CDI-) with random initial conditions are plotted. 



region of Figure lA), in which the outcome of the com- 
petition depends on the initial conditions (Figure IE). 
In this parameter regime, both of the extinction states 
are linearly stable. A plot of the Lyapunov function [36] 
(Figure IB) indeed shows two minima for the phase space 
of the system, corresponding to possible stable extinction 
states. 

The above results suggest that the tradeoff between 
inhibition strength, due to the production of toxins by 
CDI, and relative growth advantage, due to the metabolic 
load associated with CDI, determines the possible out- 
comes of a two-species competition. This conclusion is 
in agreement with a previous experimental study con- 
cerning bacteriocin production [49], where initial con- 
centrations of a competing population determined the 
extinct species. Although a diffusible toxin instead of CDI 
was employed in the experiment, it clearly supports our 
modeling results. In another experimental report, it was 
shown that a CDI+ strain is indeed capable of driving a 
CDI- strain toward extinction in a well-mixed environ- 
ment [14]. In the future, it will be valuable to design 



experiments to systematically investigate the impacts of 
initial concentrations and metabolic load on the outcome 
of CDI-based species competition. 

The dependence of competition outcomes on initial 
conditions in certain parameter regions may have addi- 
tional important implications: It indicates that spatial 
aggregation and localized patterns may occur in space 
with appropriate initial conditions when the diffusion of 
the system is small. 

Phase diagram for the existence of spatial patterns 

In a preliminary search for the existence of stable patterns 
in one spatial dimension, we performed simulations with 
random initial conditions for a specific set of the growth 
advantage {a/p), inhibition (ci = c/3), and diffusion con- 
stants {Du = Dy = D), Figure 2 illustrates three time 
courses in which the concentration distribution of species 
u shows stable localized patterns (spatial distribution of 
species v is anti-correlated with that of u and presented 
in Additional file 1: Figure Sl-3). This phenomenon orig- 
inates from the stability of both extinction states of our 
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Figure 2 Representative time course evolutions of the 
competing two-species population in one dimensional space 
with a growth advantage {oc/p) of 3.5 and an inhibition {c^ ) of 

2.1. Coexistence of the CDI+ and CDI- species may appear in certain 
parameter sets (Panel A and B). However, localized patterns will 
merge and eventually disappear with the increase of the diffusion 
constant for both species (Panel A:D = 0.001 , Panel B: D = 0.01 , and 
Panel C: D = 0.1 ). The same set of initial conditions are used for all of 
the panels. Here, only the spatiotemporal distributions of the 
concentration of species u (CDI-) are shown in the panels due to the 
mutual exclusion features of the two competing species. Those of the 
other species (CDI+, v) are available in Additional file 1 : Figure SI -3. 



system in a certain parameter regime. Figure 2 further 
shows that, as the diffusion constant changes, the basin 
of attraction for the states may vary: as the diffusion con- 
stant is changed from 10~^ to 10~^ and to 10~^, the small 
stripes in Figure 2A become unstable in Figure 2B and all 
stripes become unstable in Figure 2C. Therefore, diffusion 
tends to destroy possible species aggregations within the 
population. 

To systematically explore possible outcomes of the com- 
petition in space and to understand the conditions under 
which corresponding patterns emerge, we computation- 
ally examined possible steady states of the model as we 
varied the relative growth advantage (a/p), inhibition 
strength (ci), and species motility (D) for different ini- 
tial conditions. In principle, a parameter search can be 
implemented by repeatedly simulating our mathematical 
model (Eq. 2) with massive random initial conditions for 
each possible parameter combination and then by ana- 
lyzing the corresponding steady states of the system. To 
reduce the computational cost of the problem, we chose 
a representative set of initial conditions for simulation in 
practice. Here, the initial condition sets consisted of two 
non-overlapping domains with each solely occupied by 
one of the two species at the systems carrying capacity, 
i.e., grids 1-5 are solely occupied by u while the remain- 
ing space, grids 6-32, are fully occupied by v for the 
32-grid space (detailed in Methods and Additional file 1: 
Figure S8). 

Figure 3A shows the phase diagram for the competition 
outcome of the population obtained from our numeri- 
cal study. As illustrated in the figure, the beak regions 
outlined by the color lines are the coexistence parame- 
ter space where spatial localized patterns may emerge and 
be stable. The colors of the lines (green, cyan, and navy) 
correspond to the species' diffusion constant of 0.0001, 
0.001, and 0.01 respectively. Outside the beak regions, the 
competition of the population always lead to the extinc- 
tion of one of the two species. Although obtained using 
the representative sets of initial conditions, the phase 
boundaries were tested and confirmed by employing 10"^ 
random initial conditions for D = 0.001 with a/p g[ 3, 4]. 
To further illustrate the diversity of possible competi- 
tion outcomes, representative patterns with the highest 
likelihood from random initial conditions are shown in 
Figure 3B-G, where the inhibition strength is varied 
across the phase diagram (ci = 1.7,1.8,1.9,2.0,2.1,2.3) 
for a fixed growth advantage {a/P = 3.5) and diffusion 
constant {D = 0.001). 

It is important to notice that the spatial coexistence 
region is within the parameter region where both extinc- 
tion states remain stable (shaded region), indicating that 
the emergence of patterns is not due to a diffusion driven 
instability of the homogeneous states like Turing patterns 
[35,50]. In addition, as every spatial aggregate of the two 
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Figure 3 Phase diagram and representative patterns of the competing population. (A) Phase diagram. The grey shaded region corresponds 
to the case when both CDI+ or - extinctions are possible when the population is well mixed, identical to the grey region in Figure 1 A. The green, 
blue and navy lines are the phase boundaries between spatial coexistence (beak regions) and extinction (outside of beak regions) when the rate of 
population diffusion is 0.0001, 0.001, and 0.01 respectively. (B-G)The most likely localized patterns of the competing population with different 
growth advantages and inhibitions. 1 0^ random initial conditions are used here as we move across the phase diagram with a growth advantage of 
3.5, a diffusion constant of 1 0~^, and an inhibition varying from 1 .7 to 2.3. (B) For an inhibition of 1 .7, patterns may be present depending on initial 
conditions, but the most likely state is the extinction of species v. (C-D) For an inhibition of 1 .8 and 1 .9, the most likely pattern changes to a single 
stripe configuration. (E-F) For an inhibition of 2.0 and 2.1 , the most likely pattern changes to a double stripe configuration. (G) For an inhibition of 
2.3, only the extinction states are observed, with the most likely state being the extinction of u. 



species within a stable pattern is locally similar to the 
extinction states of the system in the well-mixed case, the 
stability of both extinction states appears to be a neces- 
sary condition for the existence of spatial patterns. Thus, 
the tradeoff between relative growth advantage and inhi- 
bition strength continues to play a key role in determining 
possible outcomes of competition in space. 

Statistics of localized patterns 

Through our computational survey of competition out- 
comes, we have found that diverse patterns may emerge 
from the competing population for a single parameter set: 
As illustrated in Figure 4, the competition of the popula- 
tion can possibly produce single-stripe, double-stripe, and 
multi-stripe localized patterns as well as the extinction 
of either species for a growth advantage of 3.5, inhibi- 
tion constant of 2.2, and a diffusion constant of 10~^. 
Although various patterns were observed, their preva- 
lence was different: In this specific case, double, single and 
multi-stripe patterns constitute the vast majority of pos- 
sible outcomes with a relative abundance of 51.5%, 23.7%, 
and 23.3% respectively. 

To reveal how key system parameters influence the 
occurrences of possible patterns, we performed the simu- 
lations of our model with 10"^ random initial conditions for 
a varying inhibition strength (relative growth advantage 
and diffusion constant are fixed) and analyzed the corre- 
sponding statistics. As shown in Figure 5A, competition 
of the population always gives rise to extinction (species 
v) when the inhibition strength is weak (ci = 1.4), corre- 
sponding to the left of the coexistence region in Figure 3A. 



However, as the inhibition strength increases, the system 
enters the coexistence parameter space and various pat- 
terns are observed with different relative occurrences. The 
competition outcome of the population becomes exclu- 
sively extinction (mostly of species u) when the inhibition 
is beyond a threshold (the right boundary of the coexis- 
tence region). Accordingly, the most probable pattern is 
changed from extinction to single-stripe patterns and later 
to double-stripe patterns before eventually returning back 
to extinction (Figure 3B-G). The observed alteration of 
pattern diversity and relative abundance is due to the fact 
that the orchestration of the tradeoff between inhibition 
and growth advantage is mandatory for the spatial coex- 
istence of the two species and imbalance of the tradeoff 
results in the increase of species extinction. In addition to 
the relative abundance, the widths of stripes of localized 
patterns also change with the inhibition strength as shown 
in Table 1 (and Additional file 1: Table S3). In contrast 
to the relative abundance of pattern occurrences, pattern 
widths are less subjective to parameter changes, mainly 
due to the requirement of minimal size of localized pat- 
terns for stabilization and the mutual exclusion feature of 
the two species. 

Similar to the inhibition strength, the relative growth 
advantage can also considerably alter the localized pat- 
terns of the competing population. As illustrated in 
Figure 5B, the relative abundances of stable patterns vary 
with the growth advantage with the most probable pattern 
changing from extinction to double-stripe to single-stripe 
and finally to extinction again. This trend of pattern alter- 
ations is opposite to the case of varying inhibition strength 
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Figure 4 Occurrence of localized patterns of the competing population with a growth advantage of 3.5, an inhibition of 2.2 and a 
diffusion constant of 10~'^. Various patterns may emerge in tiiis parameter set, among wliicli double-stripe patterns are tine most dominant 
(51 .5%), followed by single-stripe (23.7%) and then multi-stripe (23.3%) patterns. Only about 1 .5% of initial conditions result in extinction for one of 
the species. Here, 1 0^ runs with random initial conditions are used to generate the statistics. 



(Figure 5A), This is primarily due to the opposite impacts 
of growth advantage and inhibition in shaping population 
structures. The statistics of the width of stripes of stable 
patterns with different growth advantages are shown in 
Table 2 (and Additional file 1: Table S4). 

The statistics of the stripe patterns further show that 
the GDI- strain tends to have larger cluster sizes than 
the CDI+ strain on average within the coexistence region. 
In fact, for all coexistence statistics computed, the GDI- 
strain occupies a greater percentage of the space than the 
GDI+ strain. The disparity of the stripe widths for GDI+ 
and GDI- is primarily due to the non-local nearest neigh- 
bor interactions: To have a sustained domain, GDI- cells 
must build a buffer region to protect themselves when 
they are surrounded by GDI+ cells; in contrast, when 
GDI+ cells are surrounded by GDI- cells, no buffer lay- 
ers are needed. Therefore, there is an increase in colony 
width for GDI- colonies, compared to GDI+ colonies. 
Importantly, the stripe width disparity does not exist in 
systems with only on-site interactions as employed previ- 
ously in many competition studies [29-34]. Instead, both 
species will occupy the space equally on average. From an 



ecological perspective, the disparity of localized patterns 
shows the importance of non-local interactions in deter- 
mining the species aggregation of bacterial communities 
and may offer new insights into the diversity of microbes. 

Perturbative expansion for patterned states 

To gain insight into the exact requirements for inhibi- 
tion and growth advantage necessary to produce spatial 
coexistence of the population, we have employed a pertur- 
bation approach in addition to the numerical results pre- 
sented above. The perturbation approach is built around 
finding a stable spatially inhomogeneous steady state 
for the model with D = 0 and then considering how 
the state is perturbed when slow diffusion is allowed 
(see Additional file 1: Section 2 for details). 

When D = 0, each equation in our model (Eq. 2) at 
steady state becomes the product of two linear equations. 
Thus, there are 2^^ possible steady states in total, where 
N is the number of grid points. To determine the pos- 
sibility of coexistence, we consider the simplest possible 
steady state with coexistence and aggregation, in which 
species u and v each form a single domain at the respective 
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Figure 5 Statistics of localized patterns as a function of growth 
advantage and inhibition. (A) Relative abundance of spatial 
patterns as a function of CDI inhibition (ci) for a given growth 
advantage {a/^ = 3.5). The most likely steady pattern transits from 
CDI+ extinction to single to double stripe, with a increasingly high 
prevalence of multi-stripes, as CDI inhibition increases. It makes a 
sharp change back to extinction (CDI-) as the parameters exit the 
spatial coexistence region. (B) Relative abundance of spatial patterns 
as a function of growth advantage {a/^) for a given CDI inhibition 
(ci = 2.0). An opposite trend of pattern statistics is observed 
compared with panel A. 



carrying capacities with two single-grid-point transition 
layers connecting the two domains in space. By plugging 
this trial state into our model, we have: 



Ui = a 

Ui = a 

Ui = 0 
Ui = a 



3ci Vi = 0 / 



N 
~2 



N 
N 



V,- = 0 i = N 



.N-l 



Table 1 Summary of statistics for different inhibition 
values (ci ) using a growth advantage {oc/P) of 3.5 and a 
diffusion constant (D) of 0.001 



^1 


Olliyic all iptr 

width 


width 


IVIUIll~3iripc 

width 


1 .4000 


N/A 


N/A 


N/A 


1.5000 


8.9207±0.2727 


3.9583±1.2969 


0.0000 


1 .6000 


9.1718±0.3142 


4.1944±2.1219 


2.5625±1.3797 


1.7000 


9.0973±0.4039 


4.1153±2.0673 


2.4557±1.4078 


1.8000 


8.9883±0.5083 


4.0300±2.1113 


2.3128±1.3519 


1.9000 


8.8421 ±0.6373 


3.9168±2.0163 


2.2297±1.2819 


2.0000 


8.5869±0.8970 


3.741 9±1. 9762 


2.0780±1.2299 


2.1000 


8.1 909±1. 2791 


3.4835±1.8674 


1.9343±1.1349 


2.2000 


6.7693±2.1881 


3.0211 ±1.5841 


1.8590±0.9125 


2.3000 


N/A 


N/A 


N/A 



Widths are the average size of domains for species u in a system of length 
{L=^0) with the standard deviations shown. 



for which we have assumed that two transition layers 
occur at / = N/2 and / = N for simplicity (our analysis 
can be similarly applied for the transition layers being any- 
where in the space). Note that Eq. 4 above corresponds to 
Eqs. S9-10 and S12 in Additional file 1. Through linear sta- 
bility analysis, we found that this proposed state is indeed 



Table 2 Summary of statistics for different growth 
advantage values {oc/P) using an inhibition (ci ) of 2.0 and 
a diffusion constant (D) of 0.001 





Single stripe 
width 


Double stripe 
width 


Multi-stripe 
width 


3.2000 


N/A 


N/A 


N/A 


3.3000 


6.6860±2.2685 


3.0043±1.5822 


1.8632±0.9081 


3.4000 


8.2396±1.2199 


3.4892±1.8759 


1.9354±1.1538 


3.5000 


8.6220±0.8392 


3.7115±1.9583 


2.0711 ±1.2085 


3.6000 


8.8299±0.6749 


3.8809±2.0082 


2.1 776±1. 2621 


3.7000 


8.941 5±0.5602 


3.9947±2.0408 


2.2990±1.3339 


3.8000 


9.0544±0.4490 


4.0685±2.0718 


2.3592±1.3373 


3.9000 


9.1103±0.3828 


4.1275±2.0864 


2.4225±1.4374 


4.0000 


9.1619±0.3433 


4.1917±2.0968 


2.5129±1.3891 


4.1000 


9.1953±0.2962 


4.1774±2.1287 


2.5104±1.4332 


4.2000 


9.2119±0.2598 


4.2721 ±2.1 072 


2.5651 ±1.2789 


4.3000 


9.231 0±0.2421 


4.1 741 ±2.3164 


0.0000 


4.4000 


8.9557±0.2050 


0.0000 


0.0000 


4.5000 


8.9623±0.2083 


0.0000 


0.0000 


4.6000 


N/A 


N/A 


N/A 



(4) 



Widths are the average size of domains for species u in a system of length 
{L=^0) with the standard deviations shown. 
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linearly stable in the case ofD = 0 (see Additional file 1: 
Section 2). 

When the system has a nonzero diffusion {D ^ 0), 
we would like to determine if some perturbed form of 
the above steady state remains stable. This can be imple- 
mented by performing a perturbation expansion in A i.e., 
Ui ~ uoi + Duli + D'^U2i + 0(D^) and v/ ~ vq/ + Dvu + 
^^'^2i + 0(D^), where uoi and voi are the steady state val- 
ues presented in Eq. 4 {ui and v^). To have a stable state, 
it is mandatory to have all Ui and v/ remain positive and 
bounded by the respective carrying capacities, which gives 
rise to a necessary condition for stability with £) 7^ 0 as 
(see Additional file 1: Section 2) 

1 + ci < ^ < 1 + 2ci (5) 

P 

which corresponds to Eq. S20 in Additional file 1. The 
necessary condition also allows us to estimate the coexis- 
tence phase boundary in the limit as £) ^ 0, as shown in 
Additional file 1: Figure S4 

Although we only considered the perturbation results 
for two one-grid-point transition layers, our results can 
be generalized to larger transition layers. The simplicity 
of the one grid point case, however, provides an ana- 
lytical bound between growth advantage and inhibition 
strength for the existence of stable localized patterns. As 
the diffusion constant of the model approaches to zero, 
we found that our analytical calculation from Eq. 5 indeed 
represents a good approximation of the phase boundaries 
obtained from numerical simulations (see in Additional 
file 1: Figure S4). 

Two dimensional patterns 

Although we have primarily focused on the exploration 
of spatial population structures in one dimension thus 
far, coexistence and localized patterns can be present in a 
high dimensional space as well. To demonstrate our idea, 
we extended our model (Eq. 2) for the two dimensional 
case where both species diffusion and nearest neighbor 
interactions are over a square grid. Similar to the one 
dimensional case, an inhibition constant C2 = c/5 was 
introduced to reflect the five locations for possible GDI 
(one on-site and four nearest neighbors). 

To illustrate the plausibility of competition-induced spa- 
tial patterns, we examined possible outcomes of the popu- 
lation with the alteration of the strength of GDI inhibition 
while keeping the relative growth rate and diffusion con- 
stant fixed. Accordingly, two sets of representative initial 
conditions are employed: one set of initial conditions 
having a cluster of u that is centered in the space and 
surrounded by v, and the other set having a cluster of 
V centered in the space and surrounded by u (detailed 
in Methods and Additional file 1: Figure S9). Figure 6A 
shows the counts of stable coexistence patterns for both 



sets of the initial conditions with the green inverted tri- 
angles corresponding to the first set and the orange tri- 
angles corresponding to the second. It is clear that the 
counts of stable patterns are different for the two differ- 
ent initial condition sets even with the same parameter 
set, which is due to the asymmetry of the interactions 
between the GDI+ and GDI- species. By summing up 
the two counts, we found that the overlapping region of 
Figure 6A (blue region) corresponds to the parameter 
space where the population competition has the largest 
probability of forming coexisting patterns, corresponding 
to the least chance of species extinction, for the sets of 
initial conditions tested. 

With the exploration of parameter space, we further 
illustrated spatially localized patterns in two dimensional 
space in Figure 6B-D where representative spatial dis- 
tributions of species u are shown (see Additional file 1: 
Figure S5-7 for species v). The average amount of space 
occupied by species u decreases with the increase of the 
inhibition. As in the one dimensional case, the space occu- 
pied by u is larger than that of v in the coexistence region 
for the parameters explored. Our results demonstrate that 
the two dimensional model is indeed capable of producing 
coexistence for the two species with a rich set of possible 
patterned steady states. 

Conclusions and discussion 

With a two-species population model, we have computa- 
tionally investigated the dynamics and competition out- 
comes of a bacterial population with contact-dependent 
inhibition in different settings. We found that the tradeoff 
between the strength of inhibition via direct cell contact 
and relative growth advantage associated with metabolic 
burden are of central importance in determining the out- 
come of bacterial competition. In the well-mixed case, 
two competing species are mutually exclusive and their 
competition always leads to the extinction of one of the 
two species depending on the inhibition-growth tradeoff 
as well as initial conditions. In contrast, coexistence and 
localized patterns may also emerge from the competition 
of the population in the spatial case, in addition to exclu- 
sive extinctions. In addition, a statistical picture of a pop- 
ulation with GDI-based competitions has been revealed, 
including the diversity, relative abundance, and pattern 
characteristics of all possible competition outcomes. 

This study has hence expanded the spectrum of the 
current knowledge about contact dependent inhibition 
and provided a systematic view of GDI-based competition 
in bacterial populations by exploring possible outcomes 
in different settings. It also offers a quantitative calibra- 
tion on the requirements for the emergence of various 
outcomes by examining the effects of inhibition, rela- 
tive growth rate, and population diffusion. In addition, 
the population competition dynamics illustrated here 
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Figure 6 Localized patterns in two dimensional space. (A) Counts of stable localized patterns for two sets of initial conditions (see Additional 
file 1 : Figure S9). The green and orange curves (triangles) correspond to the counts for the initial conditions when u is surrounded by v and v is 
surrounded by u respectively. The overlap region (blue) suggests a parameter regime where species coexistence is stable for both sets of initial 
conditions. (B-D) Localized patterns in two dimensional space with a growth advantage of 3.5, a diffusion constant (D) of 10~^, and an inhibition of 
1 .0, 1 .05, and 1 .1 respectively. Each pattern was chosen to represent the average space occupied by CDI- for 1 00 runs with random initial conditions. 
The average percentage of space occupied by CDI- for the three inhibition values is 95%, 84%, and 78% respectively. 



increase our understanding of the complexity of bacterial 
social interactions. Furthermore, this work sheds light on 
new experimental design and tests by providing the pre- 
dictions on spatial coexistence and locaUzed patterns of 
the population structure as well as their initial condition 
dependence. 

It is worthy to note that our mathematical model 
adopted a deterministic description of population com- 
petition. However, population dynamics and competition 
outcome may differ or, in some cases, deviate dramatically 
from those of a deterministic model when we take into 
account the stochastic nature of cellular growth [51]. For 
instance, diversity of emergent patterns and their charac- 
teristics may be altered significantly when the system is 
close to its phase boundary where multiple distinct spatial 
patterns occur. In addition, contact-dependent inhibition 
may have altered modes of action in different organisms 
as the expression of the CDI genes can be constitu- 
tive or stationary phase dependent [14,25,52]. For future 
study, it will therefore be valuable to consider the growth 
phase dependence of the activities of contact dependent 
inhibitions. 



Methods 

Numerical integration 

To explore the competition outcomes of a two-species 
bacterial population, we numerically integrated our math- 
ematical model consisting of a set of coupled ordinary 
differential equations (e.g., Eqs. 2 and 3). Here, integra- 
tion was carried out using a Runge-Kutta adaptive step 
size method of order 4-5 [53,54]. The maximal error 
per step was set to 10~^ for all runs unless otherwise 
noted, and a maximal step size of 10~^ was imposed. 
To avoid negative values for any of the species con- 
centrations during integration, we decreased the time 
step until all negative values were within the maximal 
error threshold from zero and subsequently set the neg- 
ative values to zero before the next time step was taken. 
Our results for numerical integration were validated with 
random sample runs using a different maximal error 
and maximal step size as detailed in Additional file 1: 
Section 3. 

Discretization of the reaction diffusion equations for 
our system was accomplished using a Taylor expansion 
approximation to the diffusion operator. Namely: 
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cfiu u(x -\-8) — 2u(x) + u(x — 8) 



which converts the coupled partial differential equations 
of our system to a set of coupled ordinary differential 
equations. The number of coupled ordinary differential 
equations was determined by the number of grid points 
employed for the system, which for the one dimensional 
case is 32 unless otherwise noted. 

Spatial patterns 

To determine the possibility of population coexistence 
and aggregation in one dimensional space, we used initial 
conditions that involve two domains with each consisting 
solely of one of the two species at its individual carry- 
ing capacity as shown in Additional file 1: Figure S8. The 
domain size was varied between all possible sizes for each 
species with the entire space filled, which gives rise to 
31 different initial conditions used with initial condition ; 
given by: 

Ui = 0 Vi = P / = 

Ui = a Vi = 0 / = 7 + 1, . . . , 32 (7) 

where / varied from 1 to 31. This set of initial condi- 
tions was used to approximate the phase boundary for 
the coexistence of the two species. The boundary was 
further validated through the use of 10"^ random ini- 
tial conditions near the boundary for a/^ ^[3,4] and 
D = 10-^ 

The statistics for patterns within the coexistence region 
were collected through the simulation of 10"^ random 
initial conditions until convergence was reached. We 
considered the system to have converged if the maxi- 
mal difference for all grids averaged below 10~^ times 
the largest carrying capacity in the system for 10^ 
time steps. The stripe widths were calculated by deter- 
mining the dominant species at each grid point. Pat- 
terns were then classified according to the number of 
distinct stripes and widths present in the converged 
state. 

We further extended our model to give illustrative 
examples of possible pattern formation in two dimen- 
sional space. To explore parameter space, we utilized a set 
of initial conditions in which each species formed a spa- 
tial aggregate at the center of space surrounded by a sea 
of its competitor. For example, a square of 4 grid points 
of species u at its carrying capacity was surrounded by 
species v at its carrying capacity. The interior species size 
ranged from 4 grid points to 256 grid points with the total 
system size set at 1024 grid points (see Additional file 1: 
Figure S9). 



Parameter space 

For the well-mixed case, the phase diagram (Figure lA) 
was determined using linear stability analysis. The Lya- 
punov function [36] in the either extinction region was 
constructed using a = 2, = 1, and c = 2 (Figure IB). 
The representative time course for the parameter regions 
employed a = 2, ^ = 1, and c = 0.5 (Figure IC), a = 0.5, 
y6 = 1, and c = 2 (Figure ID), and a = 2, y0 = 1, and c = 2 
(Figure IE). 

In Figure 2, an identical set of initial conditions was 
used to produce the spatial patterns. The growth advan- 
tage (a/p) and inhibition (ci) were held fixed at 3.5 and 2.1 
respectively, while D was varied between 10~^, 10~^, and 
10-1 for Fi gure 2A-C respectively. 

To systematically explore the tradeoff between rela- 
tive growth advantage and inhibition in the small dif- 
fusion limit, we searched over parameter space in the 
one dimensional model with a/fi g[ 1, 5], ci g[ 0, 4], and 
D e {lO""^, 10"^, 10~^} to generate the phase diagram in 
Figure 3 A. The growth advantage and inhibition strength 
were sampled every 10"^ here. For Figure 3B-G, a growth 
advantage {a/fi) of 3.5 was used with a varying inhibition 
strength (ci), given by 1.7, 1.8, 1.9, 2.0, 2.2, 2.3 respectively 
to find the most likely steady state out of 10^ random 
initial conditions. 

In Figure 4, we used a growth advantage {a/fi) of 3.5, an 
inhibition constant (ci) of 2.2, and a diffusion constant of 
10~^. The percentages were taken from 10"^ random ini- 
tial conditions in which the steady states were classified 
according to single stripe, double stripe, multi-stripe, or 
extinction. 

As a summary of the types of patterns observed across 
the coexistence region in the parameter space. Figure 5 
shows the relative abundance of patterns when varying 
growth advantage and inhibition. In Figure 5A, the growth 
advantage and diffusion constant are fixed at 3.5 and 0.001 
respectively while the inhibition is varied between 1.4 and 
2.3 with 0.1 increments. In Figure 5B, the inhibition and 
diffusion constant are fixed at 2.0 and 0.001 respectively 
while the growth advantage is varied between 3.2 and 4.6 
with 0.1 increments. 

For the two dimensional model, a/fi and D were kept 
fixed at 3.5 and 10~^ respectively while C2 was sampled 
every 10-1 within [0, 3]. Fi gure 6A was generated using a 
set of initial conditions with each species surrounded by a 
sea of its competitor at carrying capacity (see Additional 
file 1: Figure S9). Figure 6B-D are possible steady state 
patterns using an inhibition constant (02) of 1.0, 1.05, and 
1.1 respectively. 

Analytical results 

Analytical results for the phase boundaries were deter- 
mined using asymptotic expansions and linear stability 
analysis (see Additional file 1: Section 2 for details). For an 
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introduction to asymptotic expansions and linear stability, 
see [55]. The existence of a stable coexistence state for 
slow diffusion in discrete systems was proposed in [29] for 
a two grid point system through the use of a perturbation 
theorem that remains valid for any number of grid points. 
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Additional file 1: Supplementary Information. 
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